rm(list = ls())
warnings()
library(fOptions)

logging = 1 #This variable is to give the user an option to print out checking parameters
#such as v_flag and v_fi, 1 = printing, 0= not printing

########################################################
#I. MINISUM FUNCTION----THEOREM 1---INDEPENDENT PROGRAM#
########################################################
f_h<-function(sigma,P){
  #Function f_h computes the h function: h(sigma)=C_BS(X1, sigma)-C_BS(X2,sigma) 
  #required for the theorem 1 and function f_minisum
  Call_mkt_1=unlist(P[1])
  Call_mkt_2=unlist(P[2])
  maturity_dt=unlist(P[3])
  X1=unlist(P[4])
  X2=unlist(P[5])
  S=unlist(P[6])
  risk_free=unlist(P[7])
  dividend=unlist(P[8]) #cost of carry  or dividends, take them zero in most cases
  #in the function help it is explained that to get BS formula we have to set b=risk_free
  #to see the content of the function enter GBSOption
  
  Call_BS_1=GBSOption(TypeFlag = "c", S = S, X = X1, Time =maturity_dt, r = risk_free,  
  b =risk_free, sigma = sigma)
  Call_BS_2=GBSOption(TypeFlag = "c", S = S, X = X2, Time =maturity_dt, r = risk_free,  
  b =risk_free, sigma = sigma)

  h_out=Call_BS_1@price-Call_BS_2@price
 
  if (logging==1) { 
    #to print h_out if logging=true
    cat('h_out=',h_out,'\n')
  }

  return(h_out)
}

##########################################
f_f<-function(sigma, P){
  #f_f function computes function f(sigma) defined in the paper, its minimum gives minisum
  Call_mkt_1=unlist(P[1])
  Call_mkt_2=unlist(P[2])
  maturity_dt=unlist(P[3])
  X1=unlist(P[4])
  X2=unlist(P[5])
  S=unlist(P[6])
  risk_free=unlist(P[7])
  dividend=unlist(P[8]) # cost of carry or dividends,take them zero in most cases
  # in the function help it is explained that to get BS formula we have to set b=risk_free
  # to see the content of the function enter GBSOption

  Call_BS_1=GBSOption(TypeFlag = "c", S = S, X = X1, Time =maturity_dt, r = risk_free,  
  b =risk_free, sigma = sigma)
  Call_BS_2=GBSOption(TypeFlag = "c", S = S, X = X2, Time =maturity_dt, r = risk_free,  
  b =risk_free, sigma = sigma)

  f=abs(Call_mkt_1 - Call_BS_1@price)+abs(Call_mkt_2-Call_BS_2@price)
  return(f)
} 

##########################################
#NEED TO EDIT TO INCLUDE CURRENT DATE IN PARAMETER INPUT AND OUTPUT
f_minisum<-function(ls_ParamOptions ){
  # computes both the minisum function , i.e. min of the f_function defined above, and
  # the arguments of f_function  at which the function reaches min, i.e. sigma_star
  #requires functions: f_f and f_h
  #ls_ParamOptions=list(Call_mkt_1,Call_mkt_2, maturity_dt, X1,X2,S,risk_free, dividend) 

  Call_mkt_1=unlist(ls_ParamOptions[1])
  Call_mkt_2=unlist(ls_ParamOptions[2])
  maturity_dt=unlist(ls_ParamOptions[3])
  X1=unlist(ls_ParamOptions[4])
  X2=unlist(ls_ParamOptions[5])
  S=unlist(ls_ParamOptions[6])
  risk_free=unlist(ls_ParamOptions[7])
  dividend=unlist(ls_ParamOptions[8])

  #Current date and pair category are not required for calculating minisum but needed for output
  current_date=unlist(ls_ParamOptions[9]) 
  pair_category=unlist(ls_ParamOptions[10])

  length_data=length(Call_mkt_1)

  sigma_star<-vector("numeric",length = length_data)
  flag<-vector("numeric", length = length_data)
  Minisum<-vector("numeric", length = length_data)
  Minisum_percentage<-vector("numeric", length = length_data) 
    #Minisum percentage, which equals Minisum / Stock Price
  fi<-vector("numeric", length = length_data)

  for (i in 1:length_data) {
    sigma_star[i]=0
    Minisum[i]=0
    Minisum_percentage[i]=0
    flag[i]=0 # to check that one of the conditions below was satisfied, 
    #in that case flag[i] becomes 1
 
    ls_ParamOptions_i=list(Call_mkt_1[i],Call_mkt_2[i], maturity_dt[i], X1[i],X2[i],
    S[i],risk_free[i], dividend[i]) 
  
    fi[i]=(log(X1[i]*X2[i]/(S[i]*S[i]))/maturity_dt[i])-2*risk_free[i]

    sigma_imp1=GBSVolatility(price=Call_mkt_1[i], TypeFlag="c", S=S[i], X=X1[i], 
    Time=maturity_dt[i], r=risk_free[i], b=risk_free[i], tol=0.0001, maxiter=100)
    sigma_imp2=GBSVolatility(price=Call_mkt_2[i], TypeFlag="c", S=S[i], X=X2[i], 
    Time=maturity_dt[i], r=risk_free[i], b=risk_free[i], tol=0.0001, maxiter=100)

    if (logging==1) { 
      cat('sigma_imp1=', sigma_imp1,'\n')
      cat('sigma_imp2=', sigma_imp2,'\n')
    }

    f1=f_f(sigma_imp1,ls_ParamOptions_i)
    f2=f_f(sigma_imp2,ls_ParamOptions_i)

    #Compare sigma_imp1, sigma_imp2 and fi_sqrt to determine value of minisum at sigma_star
    if(fi[i]<=0){
      Minisum[i]=f2
      sigma_star[i]=sigma_imp2
      flag[i]=1.1+flag[i]
      next # forces to start a new iteration of the loop
    } # end if(fi[i]<=0)

    if(fi[i]>0) { 
      fi_sqrt=sqrt(fi[i])
      if( (fi[i]>0) && (fi_sqrt<= min(sigma_imp1, sigma_imp2 ))){
	   Minisum[i]=f2
	   sigma_star[i]=sigma_imp2
	   flag[i]=1.2+flag[i]
	   next
      }# end if(fi[i]>0 && fi_sqrt<= min(sigma_imp1, sigma_imp2 ))

      if( (fi[i]>0) && (sigma_imp1 < fi_sqrt) && (fi_sqrt < sigma_imp2)){
        Minisum[i]=min(f1,f2)

	  #cat('Inside function f_minisum: now will be computing sigma_star=argmin( h(), h() \n')
	  #cat('sigma_imp1=',sigma_imp1  ,  'sigma_imp2=',sigma_imp2,'\n')
      
        if(f_h(sigma_imp1,ls_ParamOptions_i)<f_h(sigma_imp2,ls_ParamOptions_i)) {
          sigma_star[i]= sigma_imp1} 
          else {
            sigma_star[i]= sigma_imp2} #end if( f_h(sigma_imp1,ls_ParamOptions_i)<....

          if (logging==1) {
            cat('Inside function f_minisum' ,'\n')
          }

	    flag[i]=1.3+flag[i]
	    next
	}# end if( (fi>0) && (sigma_imp1 < fi_sqrt) && (fi_sqrt < sigma_imp2)){

	if((fi[i]>0) && (sigma_imp2 < fi_sqrt) && (fi_sqrt < sigma_imp1) ){
	  Minisum[i]=f_f(fi_sqrt,ls_ParamOptions_i)
	  sigma_star[i]=fi_sqrt
	  flag[i]=1.4+flag[i]
	  next
	}#end if((fi[i]>0) && (sigma_imp2 < fi_sqrt) && (fi_sqrt < sigma_imp1) )

      if((fi[i]>0) && (fi_sqrt >max(sigma_imp1,sigma_imp2))){
        Minisum[i]=f1
	  sigma_star[i]=sigma_imp1
	  flag[i]=1.5+flag[i] 
	  next
      }# end if(((fi[i]>0) && (fi_sqrt >max(sigma_imp1,sigma_imp2) )  )
    }

    stop('Error in the function f_minisum: this line must not be reached ')
  }# end for (i in 1:length(Call_mkt_1)

  #To calculate Minisum Percentage
  for(i in 1:length(Call_mkt_1)) {
    Minisum_percentage[i]=100*Minisum[i]/S[i]
  }
 
  ls_out=list(Call_mkt_1, Call_mkt_2, maturity_dt, X1,X2,S,risk_free, dividend, 
  Minisum,sigma_star, flag, fi,Minisum_percentage, current_date, pair_category)
  
  return(ls_out)
}# end f_minisum

###############################################
f_ReadData<-function(file_location_name){
  #This function reads data, which location and name are given in the string file_location_name
  #put data into several vectors and return a list of those vectors

  print('f_ReadData: before data reading')
  d<-read.table(file_location_name, header=TRUE,sep=",",na.string=".",as.is=TRUE)
  print('f_ReadData: after data reading')

  length_data=length(d$Call_mkt_1)
  print(length_data)

  v_Call_mkt_1<-vector("numeric",length = length_data)
  v_Call_mkt_2<-vector("numeric",length = length_data)
  v_X1<-vector("numeric",length = length_data)
  v_X2<-vector("numeric",length = length_data)
  v_maturity_dt<-vector("numeric",length = length_data)
  v_S<-vector("numeric",length = length_data)
  v_risk_free<-vector("numeric",length = length_data)
  v_dividend<-vector("numeric",length = length_data)
  v_current_date<-vector("numeric",length = length_data)
  v_pair_category<-vector("numeric",length = length_data)

  for (i in 1:length_data) { 
    #Require X1<X2
    if(d$X1[i]<=d$X2[i]){
      v_Call_mkt_1[i]=d$Call_mkt_1[i]
      v_Call_mkt_2[i]=d$Call_mkt_2[i]
      v_X1[i]=d$X1[i]
	v_X2[i]=d$X2[i]
    }  

    else { # in case X1>X2 we switch the options
      v_Call_mkt_2[i]=d$Call_mkt_1[i]
	v_Call_mkt_1[i]=d$Call_mkt_2[i]
	v_X2[i]=d$X1[i]
	v_X1[i]=d$X2[i]
    }# end if(d$X1[i]<=d$X2[i])

    v_maturity_dt[i]=d$maturity_dt[i]
    v_S[i]=d$S[i]
    v_risk_free[i]=d$risk_free[i]
    v_dividend[i]=d$dividend[i]
    v_current_date[i]=d$current_date[i]
    v_pair_category[i]=d$pair_category[i]

    ls_ParamOptions=list(v_Call_mkt_1, v_Call_mkt_2, v_maturity_dt,
    v_X1,v_X2,v_S,v_risk_free, v_dividend,v_current_date,v_pair_category) 

  } #end for (i in 1:length_data)

  if (logging==1) { 
    print('f_ReadData: ls_ParamOptions=')
    print(ls_ParamOptions)
  }

  return(ls_ParamOptions)
}#end f_ReadData<-function(d)

###############################################
f_output_minisum<-function(ls_result_minisum,s_output_file){
  v_Call_mkt_1=unlist(ls_result_minisum[1])
  v_Call_mkt_2=unlist(ls_result_minisum[2])
  v_maturity_dt=unlist(ls_result_minisum[3])
  v_X1=unlist(ls_result_minisum[4])
  v_X2=unlist(ls_result_minisum[5])
  v_S=unlist(ls_result_minisum[6])
  v_risk_free=unlist(ls_result_minisum[7])
  v_dividend=unlist(ls_result_minisum[8])
  v_Minisum=unlist(ls_result_minisum[9])
  v_sigma_star=unlist(ls_result_minisum[10])
  v_flag=unlist(ls_result_minisum[11])
  v_fi=unlist(ls_result_minisum[12])
  v_Minisum_percentage=unlist(ls_result_minisum[13])
  v_current_date=unlist(ls_result_minisum[14])
  v_pair_category=unlist(ls_result_minisum[15])
  
  #m_result_minisum=cbind(v_Call_mkt_1,v_Call_mkt_2,v_maturity_dt,v_X1,v_X2,v_S,
  #v_risk_free,v_dividend,v_Minisum,v_sigma_star,v_flag, v_fi, v_Minisum_percentage )
  
  m_result_minisum <-data.frame(cbind(v_X1,v_X2,v_Call_mkt_1,v_Call_mkt_2,v_maturity_dt,v_S,
  v_risk_free,v_dividend,v_current_date,v_pair_category,v_Minisum, v_Minisum_percentage ), 
  stringsAsFactors = default.stringsAsFactors())
 
  write.table(m_result_minisum , file = s_output_file, sep = ",", col.names = NA, qmethod = "double")
  print(m_result_minisum )

  return(0)
} #end f_output(ls_result_minisum,s_output_file )


#########################################
###I.2. MAIN PART of the PROGRAM MINISUM#
#########################################
# read data
s_input_file="C:/Users/diepnguyen/Documents/MSF- Clark/Implied Binomial Tree/Research/R/OptionPairs.csv"
ls_data=f_ReadData(s_input_file)
print(ls_data)

#compute minisum
ls_result_minisum=f_minisum(ls_data)
print(ls_result_minisum)

#print out outcome
s_output_file = "C:/Users/diepnguyen/Documents/MSF- Clark/Implied Binomial Tree/Research/R/Result_Minisum.csv"
f_output_minisum(ls_result_minisum,s_output_file )

getwd() #See directory to the output file
